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Abstract — For  more  than  a  decade,  the  U.S.  government  has 
been  developing  laser-based  sensors  for  detecting,  locating,  and 
classifying  aerosols  in  the  atmosphere  at  safe  standoff  ranges.  The 
motivation  for  this  work  is  the  need  to  discriminate  aerosols  of 
biological  origin  from  interferent  materials  such  as  smoke  and 
dnst  using  the  backscatter  from  multiple  wavelengths  in  the  long 
wave  Infrared  (LWIR)  spectral  region.  Through  previous  work, 
algorithms  have  been  developed  for  estimating  the  aerosol  spec¬ 
tral  dependence  and  concentration  range  dependence  from  these 
data.  The  range  dependence  is  required  for  locating  and  tracking 
the  aerosol  plumes,  and  the  backscatter  spectral  dependence  is 
used  for  discrimination  by  a  support  vector  machine  classifier. 
Snbstantial  progress  has  been  made  in  these  algorithms  for  the 
case  of  a  single  aerosol  present  in  the  lidar  line-of-sight  (LOS). 
Often,  however,  mixtures  of  aerosols  are  present  along  the  same 
LOS  overlapped  in  range  and  time.  Analysis  of  these  mixtures 
of  aerosols  presents  a  difflcnlt  inverse  problem  that  cannot  be 
successfully  treated  by  the  methods  used  for  single  aerosols. 
Fortunately,  recent  advances  have  been  made  in  the  analy¬ 
sis  of  inverse  problems  nsing  shrinkage-based  Li -regularization 
techniques.  Of  the  several  Li -regularization  methods  cnrrently 
known,  the  split  Bregman  algorithm  is  straightforward  to  imple¬ 
ment,  converges  rapidly,  and  is  applicable  to  a  broad  range  of 
inverse  problems  including  our  aerosol  unmlxlng.  In  this  paper, 
we  show  how  the  split  Bregman  algorithm  can  successfully  resolve 
LWIR  lidar  data  containing  mixtures  of  bioaerosol  simulants 
and  interferents  into  their  separate  components.  The  individual 
components  then  can  be  classified  as  bio-  or  nonbioaerosof  by  onr 
SVM  classifier.  We  iliustrate  the  approach  throngh  data  coliected 
in  field  tests  over  the  past  several  years  using  the  U.S.  Army  FAL 
sensor  in  testing  at  Dugway  Proving  Gronnd,  UT. 

Index  Terms — Aerosols,  lidar,  Li -regularization,  parameter 
estimation,  remote  sensing. 

I.  Introduction 

This  work  describes  a  generalization  of  previous 
research  [I]  on  developing  efficient  algorithms  for 
estimating  the  parameters  of  optically  thin  aerosols  in  the 
atmosphere  using  data  from  rapidly  tuned  multiple-wavelength 
long  wave  infrared  (LWIR)  lidar  [2].  The  motivation  for  this 
work  remains  the  same:  the  need  for  detecting,  locating,  track- 
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ing,  and  discriminating  atmospheric  aerosols  at  safe  standoff 
ranges  using  time- series  data  collected  at  a  discrete  set  of  CO2 
laser  wavelengths.  Our  earlier  work  considered  a  single  aerosol 
material,  producing  a  backscatter  enhancement  over  that  from 
the  natural  atmosphere.  The  goals  were  to  detect  and  track 
the  aerosol  by  means  of  estimates  of  the  range  dependence 
of  the  concentration,  and  to  discriminate  potentially  harmful 
aerosols — particularly  those  of  biological  origin — from  inter¬ 
ferents  such  as  smoke  and  dust  by  means  of  the  spectral 
dependence  of  the  backscatter.  State-of-the-art  support  vector 
machine  (SVM)  classifiers  have  been  developed  for  this  dis¬ 
crimination. 

Often,  however,  mixtures  of  aerosols  are  present  along  the 
same  lidar  line-of-sight  (LOS),  overlapped  in  range  and  time. 
For  example,  it  is  quite  possible  to  have  a  biological  or  chemical 
agent  release  from  a  munition  accompanied  by  dust  and  other 
byproducts  of  the  explosive  release.  Those  extra  materials 
can  distort  the  spectral  dependence  of  the  threat  material  and 
degrade  the  ability  of  the  sensor  to  correctly  identify  the  threat. 
The  analysis  of  aerosol  mixtures  presents  a  difficult  inverse 
problem  that  cannot  be  successfully  treated  by  the  methods 
used  previously  for  single  materials. 

Fortunately,  recent  advances  have  been  made  in  the  analysis 
of  inverse  problems  using  shrinkage-based  Li -regularization 
techniques.  We  focus  here  on  the  split  Bregman  algorithm  [3] 
because  it  is  straightforward  to  implement,  converges  rapidly, 
and  is  applicable  to  a  broad  range  of  inverse  problems  including 
our  aerosol  unmixing.  In  this  paper,  we  show  how  the  split 
Bregman  algorithm  can  successfully  resolve  LWIR  lidar  data 
containing  mixtures  of  bioaerosol  simulants  and  interferents 
into  their  separate  components.  The  individual  components 
then  can  be  classified  as  bio-  or  nonbioaerosol  by  our  SVM 
classifier. 

The  use  of  Bregman  iteration  to  solve  difficult  Li-estimation 
problems  efficiently  was  described  by  Yin  et  al.  [4].  In  that 
paper,  they  showed  the  equivalence  of  Bregman  iteration  to 
the  augmented  Lagrangian  method  for  the  minimization  of 
convex  functions  with  linear  constraints.  In  a  subsequent  paper, 
Goldstein  and  Osher  [3]  showed  that  Bregman  iteration  could 
be  efficiently  implemented  by  splitting  the  original  combination 
of  Li  and  L2  optimizations  over  a  single  variable  into  two 
separate  minimizations:  one  over  a  differentiable  quadratic 
component,  and  a  second  over  a  nondifferentiable  Li-norm 
component  with  an  objective  function  that  could  be  solved  in 
closed  form  by  shrinkage.  This  is  the  origin  of  the  term  “split 
Bregman.”  Subsequently,  [5],  [6],  it  was  recognized  that  their 
splitting  algorithm  could  be  interpreted  from  a  Lagrangian  and 
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penalty  standpoint  (augmented  Lagrangian)  analogous  to  that 
described  in  [4].  It  is  this  augmented  Lagrangian  interpretation 
of  split  Bregman  that  we  follow  here.  For  a  comprehensive 
analysis  of  the  interrelations  between  Bregman  iteration,  aug¬ 
mented  Lagrangians,  split  Bregman,  and  other  techniques,  see 
the  review  by  Esser  [7], 

The  spectral  unmixing  problem  is  of  course  well  known 
in  the  context  of  hyperspectral  (HS)  imagery,  and  hundreds 
of  papers  have  been  devoted  to  algorithms  for  performing 
it  [8].  In  the  simplest  case  when  multiple  scattering  from 
neighboring  topographical  features  can  be  ignored,  the  HS 
unmixing  problem  can  be  solved  by  least  squares  methods  if  the 
endmember  spectra  are  known.  The  latter  can  be  identified  by 
principal  components  analysis  using  training  data.  Our  problem 
is  roughly  similar  to  HS,  but  quite  distinct.  First,  whereas  HS 
images  have  typically  hundreds  of  spectral  lines  available  and 
of  the  order  of  10^  pixels  in  a  fixed  image,  our  FAL  data 
have  only  19  CO2  wavelengths  and  at  most  hundreds  of  time 
steps  and  range  cells  in  a  dynamically  changing  environment. 
Most  importantly,  we  have  no  well-defined  prior  spectra  on 
which  to  base  the  analysis;  we  must  estimate  both  the  spectra 
and  concentration  (the  analog  of  abundance  in  HS)  from  the 
same  data.  For  these  reasons,  the  lidar  unmixing  problem  is 
substantially  more  difficult  than  HS,  and  we  are  not  aware  of 
prior  publications  in  this  area. 

In  our  treatment  of  single  aerosol  estimation,  we  gave  a  rather 
detailed  discussion  of  the  lidar  measurement  and  preprocessing 
needed  to  remove  the  effects  of  the  natural  atmosphere  and  to 
deconvolve  the  long  CO2  transmitter  pulse  waveforms.  Those 
discussions  pertain  here  without  any  modification.  However,  to 
make  the  aerosol  unmixing  discussion  intelligible  to  a  broader 
readership,  in  Section  II,  we  summarize  the  main  ideas  underly¬ 
ing  rapidly  tuned  multiwavelength  backscatter  lidar  for  aerosol 
standoff  sensing.  Following  that,  we  concentrate  on  the  problem 
of  parameter  estimation  for  multiple  aerosols  by  use  of  the 
split  Bregman  method.  Accordingly,  Section  III  provides  a  brief 
summary  of  the  advantages  of  Li -regularization  over  the  more 
traditional  L2  method,  and  more  particularly,  how  this  newer, 
more  robust  method  can  be  implemented  efficiently  through 
split  Bregman.  Section  IV  describes  how  the  multi-aerosol 
estimation  is  formulated  as  a  dual  application  of  the  basic  split 
Bregman  algorithm.  In  Section  V,  we  give  some  examples  of 
processing  aerosol  mixture  data  from  both  simulated  and  actual 
release  data  collected  by  the  U.S.  Army  FAL  sensor  during 
testing  at  Dugway  Proving  Ground,  UT.  Section  VI  summarizes 
and  concludes. 

II.  Multiwavelength  Lidar  for  Aerosol  Sensing 

In  this  section,  we  review  the  main  ideas  behind  aerosol 
standoff  sensing  using  rapidly  tuned  multiwavelength  elastic 
backscatter  lidar  operating  in  the  LWIR  spectral  region.  A  more 
comprehensive  treatment  is  given  in  [1]. 

The  measurement  configuration  is  as  follows.  At  each  time 
step  a  set  of  M- wavelength-pulses  separated  in  time  by  5  ms 
is  transmitted  and  partially  backscattered  to  the  lidar  receiver. 
The  transmitter  cycles  through  the  wavelengths  in  each  set  in 
a  prescribed  pattern  chosen  either  to  emphasize  certain  spectral 


Fig.  1.  Transmitted  and  received  lidar  waveforms  at  a  single  time  step  and 
wavelength. 


features  in  the  9-11  ^m  spectral  region  accessible  to  the  lidar, 
or  to  provide  a  reasonably  uniform  wavelength  sampling  of  that 
entire  range.  Each  transmitter  pulse  waveform  is  continuously 
scattered  by  aerosols  along  the  lidar  LOS.  The  backscattered 
temporal  waveforms  are  photodetected,  electronically  hltered, 
digitized,  and  stored  at  each  wavelength  within  the  transmitted 
set  at  each  time  step. 

The  goals  of  the  processing  are  to  use  the  digitized  transmit¬ 
ted  and  received  backscatter  array  data  to  1)  decide  if  significant 
aerosol  is  present;  2)  provide  estimates  of  the  range  and  size  of 
the  aerosol  cloud;  and  3)  produce  estimates  of  the  backscatter 
spectral  dependence  for  use  by  an  aerosol  classifier.  Because 
of  the  time-series  nature  of  the  data  collection,  these  functions 
must  be  performed  sequentially  using  the  current  and  past  data 
together  with  background  data  samples  (data  not  containing 
release  aerosols)  collected  at  the  deployment  site.  These  back¬ 
ground  data  samples  are  essential  for  normalizing  the  data  that 
are  to  be  analyzed  for  the  presence  of  target  aerosols. 

Prior  to  the  main  processing  algorithms  used  to  estimate 
the  aerosol  backscatter  spectral  dependence  and  concentration 
range  dependence,  there  is  a  preprocessing  step  to  estimate  and 
remove:  1)  nonzero  baselines  in  the  transmitted  and  received 
waveforms;  2)  the  backscatter  from  the  ambient  atmosphere; 
and  3)  the  distortion  caused  by  the  transmitted  pulse  wave¬ 
forms.  The  latter  effect  is  the  result  of  a  combination  of  laser 
firing  delay,  long  CO2  pulse  duration,  and  random  shape  effects 
in  the  tail  of  the  pulse.  A  deconvolution  filter  based  on  the 
Wiener-Helstrom  method  was  used  for  the  pulse  deconvolution. 
Fig.  1  shows  these  effects  by  plotting  the  transmitted  and 
received  waveforms  versus  range  cell  from  time  step  200  at 
one  (10R38)  of  the  19  transmitted  CO2  wavelengths  of  the 
mixture  release  TD4065  from  the  U.S.  Army  FAL  sensor.  The 
natural  atmospheric  aerosol  backscatter  is  present  up  to  about 
range  cell  250.  This  structure  must  be  removed  prior  to  the 
backscatter  and  concentration  estimation.  The  smaller  structure 
around  range  cell  320  is  the  target  aerosol,  the  bacillus  BG  in 
this  case. 


WARREN  etal.-.  MULTIPLE  AEROSOL  UNMIXING  BY  THE  SPLIT  BREGMAN  ALGORITHM 


3273 


Our  model  for  the  aerosol  backscatter  at  the  lidar  receiver 
following  baseline  subtraction  and  ambient  aerosol  backscatter 
removal  is  then 

k 

Pijk  —  ^  ^  +  l  “f  ^ijk  (1) 

fe'  =  l 

where  Pijk  is  the  received  lidar  waveform  at  time  step  i,  wave¬ 
length  j,  1  <  j  <  M,  and  digitized  range  cell  k,  1  <  k  <  N, 
Tijk  is  the  corresponding  transmitted  lidar  pulse  waveform,  and 
Gijk  is  the  lidar  response  for  a  delta-pulse  transmitter,  and  Sijk 
is  an  additive  noise  term  taken  to  be  zero  mean,  independent, 
and  identically  distributed  (i.i.d.)  for  each  wavelength.  We 
model  G  as 

L 

Gijk  —  9j  ^  ^  PijlGlik  (2) 

1=1 

in  terms  of  the  aerosol  backscatter  coefficient  piji  for  material  I, 
1  <l  <  L,  Ciik  is  the  aerosol  range-dependent  concentration, 
and  gj  is  a  generally  wavelength- dependent  lidar  efficiency 
factor.  The  wavelength  dependence  reflects  that  of  the  detector 
and  optics.  In  the  absence  of  knowledge  of  g,  we  set  it  to  unity. 
The  pulse-deconvolution  filter  provides  us  with  an  estimate 
ofG 

Gijk  —  Gijk  “t”  'k^ijk  (3) 

where  riijk  is  an  additive  noise  from  the  filtered  input  noise 
e,  again  assumed  to  be  zero  mean  and  i.i.d.  The  wavelength 
dependence  of  the  backscatter  coefficients  p  is  a  composite 
result  of  the  wavelength  dependence  of  the  complex  refractive 
index  of  the  aerosol  materials  and  diffraction  effects  from  the 
particle  size  distributions  of  each  of  the  L  aerosols. 

Our  goal  in  the  data  processing  is  to  estimate  p  and  G  as 
functions  of  time,  wavelength  (p)  or  range  (C),  and  material. 
Because  of  their  appearance  in  G  as  the  sum  of  products, 
we  cannot  expect  to  obtain  estimates  having  the  correct  units 
without  further  information  and  must  settle  for  relative  esti¬ 
mates.  For  simplicity,  we  have  required  that  p  for  each  time 
and  material  be  unit  vectors  in  wavelength  lying  in  the  positive 
orthant 

M 

=  P^ji>0  (4) 

i=i 

and  Ciik  >  0. 

III.  Li -Regularization  of  Inverse  Problems 
BY  THE  Split  Bregman  Method 

It  is  well  known  that  inverse  problems  (such  as  the  retrieval 
of  signals  observed  after  some  smoothing  operation  such  as 
blurring)  are  ill-posed  in  the  sense  of  Hadamard:  small  changes 
in  the  data  due  to  noise,  for  example,  can  have  large  effects 
on  the  solution.  Otherwise  stated,  many  different  solutions  are 
compatible  with  the  same  data.  Since  most  of  these  “solutions” 
are  wildly  oscillatory,  and  therefore  not  physically  meaningful, 
some  sort  of  regularization  is  required  to  produce  useful  results. 
Historically,  L2-regularization  has  been  used  in  the  context  of 


least-squares  estimation  since  it  is  easy  to  implement  through 
ridge  regression.  Although  stabilizing  the  solution,  ridge  regres¬ 
sion  tends  to  be  unselective  in  the  features  it  suppresses  and  can 
lead  to  severely  biased  results. 

As  noted  above,  more  recent  alternatives  to  L2 -regularization 
have  replaced  the  quadratic  regularization  term  in  the  objective 
function  with  an  Li-norm  regularizing  term.  This  has  the 
effect  of  selectively  removing  features  in  the  estimate  that 
produce  the  large  fluctuations  associated  with  direct  inversion 
without  regularization.  This  selectivity  gives  Li  estimates  a 
sparse  representation  in  an  appropriate  basis  set.  The  downside 
of  Ti -regularization  has  been  the  loss  of  differentiability  of 
the  objective  function  making  implementation  more  difficult 
than  ^2- 

As  noted  above,  the  new  class  of  shrinkage-based  inver¬ 
sion  methods  such  as  split  Bregman  largely  overcomes  the 
computational  limitations  of  Li -regularization.  To  illustrate  the 
advantages  of  Ti  over  L2,  we  compare  them  on  a  simple  1-D 
simulated  example.  Our  objective  is  to  minimize  the  functionals 
Ji{u)  and  J2{u) 

ui  2  =  argmin  Ji_2(u), 

U 

Ji{u)  =  ]^\\Ku-  f\\l  +  \\\u\\i, 

J2{u)=\\\Ku-f\\l  +  ^\\u\\l  (5) 

where  Ji  2  are  the  sum  of  a  quadratic  term  quantifying  the 
fidelity  of  the  linear  model  Ku  to  the  data  /,  and  the  second 
terms  are  regularizers  using  the  Li  and  L2  norms,  respectively 

1 1/2 

.  (6) 

.  k 

We  assume  the  model  operator  K  and  the  data  to  be  known.  The 
scalar  coefficients  A  and  p  are  set  heuristically  to  balance  the 
conflicting  needs  for  fitting  the  model  and  solution  smoothness. 

For  concreteness,  we  modeled  u  as  the  sum  of  two  narrow 
Gaussian  functions  given  by 

u{x)  =  exp  [— (x  -|-  0. 1)^/0. 001]  -I-  exp  [— (cc  —  0. 1)^/0. OOl] 

(7) 

and  shown  in  Fig.  2  sampled  uniformly  over  1000  points 
between  —2  and  2.  The  data  /  were  created  by  blurring  u  with 
a  Toeplitz  matrix 

Kr^m  =  (8) 

oU 

with  r  =  0.99,  and  adding  zero-mean  independently  distributed 
(pseudo-)random  normal  noise  with  standard  error  a  =  0.004 

fn  —  ^  ^  Qn-> 

m 

g„^iV(0,iT2).  (9) 

The  simulated  data  are  shown  in  Fig.  2  together  with  the  L2 
regression  estimate  at  the  value  p  =  0.015,  and  the  Li  estimate 
computed  by  the  split  Bregman  method  defined  below  also 
using  A  =  0.015.  We  note  the  L2  estimate  is  unable  to  resolve 


iiuiii  =  y]  iwfei,  iiwii2  = 

k 
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Fig.  2.  Data  and  estimates  from  the  blurred  narrow  Gaussian  example. 

the  peaks,  whereas  the  Li  estimate  gives  a  clean  separation  of 
the  two  components. 

The  objective  function  Ji  is  nonsmooth  because  of  the  Li- 
norm  regularization  term.  This  makes  the  optimization  over  u 
difficult  by  the  usual  calculus-based  methods.  The  basic  idea  of 
split  Bregman  is  to  introduce  additional  variables  into  the  Ji 
objective  function  that  will  allow  a  convenient  separation  of  the 
minimizations  over  the  quadratic  (smooth)  data  term 

H{u)  =  ^\\Ku-f\\l  (10) 

and  the  nonsmooth  regularization  term  ||u||i.  This  is  done  in 
two  steps.  First,  the  unconstrained  problem  arg  min  H (u)  + 

U 

||u||i  is  replaced  by  the  equivalent  constrained  problem 

argminiF(w)  +  ^||(i||i  such  that  u  =  d.  (11) 

U 

Second,  the  constraint  u  =  dis  enforced  by  adding  a  Lagrange 
multiplier  term  and  an  additional  quadratic  penalty  term  to  get 
the  Lagrangian 

L{u,  d,  b)  =  H{u)  +  ^i\\d\\i  +  X{b,  u  —  d)  +  ^\\u  —  d\\l 

(12) 

where  b  is  the  Lagrange  multiplier  vector,  (a,  b)  denotes  the 
inner  product  of  vectors  a  and  b,  and  and  A  are  positive 
regularization  parameters.  In  this  formulation,  L  assumes  the 
form  of  an  augmented  Lagrangian  [4],  [5]  to  be  minimized  for 
u  and  d,  and  maximized  for  the  dual  variable  b.  In  the  Appendix, 
the  derivation  of  this  saddle  point  optimization  is  sketched.  The 
resulting  algorithm  is  given  as  Algorithm  1 .  The  function  S\  (x) 
appearing  in  Algorithm  1  is  the  shrinkage  operator  that  plays  a 
key  role  in  the  success  of  the  method  by  promoting  sparseness 
in  the  solution  through  shrinking  input  values  toward  zero.  S  is 
defined  as 

{X  —  X,x  >  X  ^ 

0,  —A  <  X  <  X  >  .  (13) 

X  +  X,x  <  —X  J 


Algorithm  1.  Split  Bregman  algorithm. 


Initialization:  d®  =  6°  =  0,  fc  =  0,  ||  Au||2  =  1 
while  ||Ait||2  >  e, 

yfc-ti  ^  xiy\K^f  +  x{d^  -  y)), 

l^k+l  =  5*:  _|_  yk+1  _  f]k+l^ 

Au  =  —  U^, 

k  ^  k  +  1, 

end 


IV.  Multi- Aerosol  Split  Bregman  Algorithm 

In  this  section,  we  generalize  the  basic  split  Bregman  method 
described  in  Section  III  to  treat  the  aerosol  unmixing  problem 
for  multiple  aerosols  within  a  single  lidar  LOS.  Our  starting 
point  is  the  M  x  N  matrix  G  =  {Gjk}  given  by  (2)  (suppress¬ 
ing  the  time-index  i)  that  represents  the  range-resolved  lidar 
data  at  a  single  time  step  after  the  preprocessing  to  remove 
the  natural  aerosol  backscatter  and  transmitter  pulse  effects  as 
summarized  in  Section  11. 

Our  goal  is  to  resolve  G  into  L  aerosol  components  using 
a  generalization  of  the  Li -regularized  split  Bregman  method. 
The  problem  then  is  to  solve 

„  min  Mp|p|i+ Ancle'll +  ;^||G-pC|||.  (14) 

||p1|  =  1,/9>0,C>0  z 


where  |a:|i  denotes  the  Li  norm  Efc|a;fe|,  and  ||||i?  indicates  the 
matrix  Frobenius  norm 


\\G-pG\\%  =  J2 


G,k-Yl  PjiCik 


1=1 


(15) 


Due  to  the  coupling  between  p  and  C  through  the  term 
II G  —  pG\W  we  do  not  have  a  convex  problem.  Our  approach 
is  to  iteratively  solve  for  p  assuming  C  and  vice  versa  through 
a  series  of  subproblems  each  of  which  is  convex  and  solvable 
by  the  split  Bregman/augmented  Lagrangian  method.  Although 
the  method  has  worked  well  on  simulated  and  actual  data,  we 
cannot  claim  optimality  for  the  overall  solution. 

To  split  the  problem  into  a  series  of  manageable  steps,  we 
introduce  the  auxiliary  variables  r  and  c  for  p  and  G,  and 
formulate  the  augmented  Lagrangian  problem 

+  +  ylk-C-f'clIL  (16) 

The  parameters  Pp,Xr  >  0  and  pc,  Ac  >  0  control  the  balance 
between  fitting  errors  and  parameter  smoothness,  and  are  set 
heuristically  using  a  combination  of  simulated  and  test  data. 
The  auxiliary  parameters  r  and  c  serve  two  functions:  they 
enforce  the  constraints  and  split  the  Li  and  L2  minimizations. 
The  parameters  b^  and  b^  are  the  Lagrange  dual  variables 
(subgradients)  that  enforce  the  constraints  p  =  r  and  C  =  c. 
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Because  the  auxiliary  variables  now  carry  the  positivity  and 
unit  norm  constraints,  the  subproblems  for  p  and  C  are  un¬ 
constrained,  and  can  be  solved  by  differentiation.  The  resulting 
updates  are 

C  =  {p^p  +  (p'^G  +  A,(c  -  6,)) 

p  =  [CC'^  +  Xriy^  {CG^  +  Xr{r  -  K))  .  (17) 

The  constrained  subproblem  for  c  becomes 

min/ic|c|i  +  ^||c- C- 6c|If  (18) 

c>0  z 

having  the  solution 

c  =  max(C'  +  bc-  /ic/Ac,0)  =  (C  +  6c  -  Mc/Ac)’''-  (19) 

The  subproblem  for  r  is  slightly  more  complicated  since  we 
require  both  ||r||  =  1  and  r  >  0.  We  have  to  solve 

min  Mpkli  +  ^lk-p-^rlll-  (20) 

||rl|  =  l,r>0  A 

The  constrained  solution  for  r  is 

_  (p  +  6c  —  Pp/Xr)^ 

||(p  +  6c-Mp/Ac)+ir  ^  ^ 

We  have  verified  that  this  solution  does  indeed  satisfy  the  KKT 
(Karush-Kuhn-Tucker)  optimality  conditions  [9].  Finally,  the 
subgradients  at  each  iteration  are  updated  as 

br  ^br  +  p  -  r, 

6c  ^  6c  +  C  -  c.  (22) 

Algorithm  2  gives  the  resulting  unmixing  algorithm.  The 
recursions  in  Algorithm  2  are  initialized  by  prior  estimates 
Pprior  of  the  backscatter,  if  available,  or  randomly.  Typically, 
10-30  iterations  are  enough  to  give  good  results. 


Algorithm  2.  Split  Bregman  algorithm  for  multi-aerosol 
unmixing. 

Initialize  p°  through  prior  estimate  Pprior  or  randomly. 

Dp  =  l,Poid  =  0 
while  Dp  >  £ 

pO  =  C°  =  6°  =  6°  =  0,  n  =  0,  (7°  =  0,  || ACHf  =  1 
while  IIACIIf  >  £ 

C'"+i  =  (p"V  +  Ac/)^\p"^G'-fAc(c'*-6^)), 
c"+i  =  (C;"+i+6^-/ic/Ac)+, 

6”+i  =  b^  +  C"+i  -  c”+\ 

AC  =  C+i  -  C", 

n  — >■  n  -f  1, 

end 

C  =  C-+MlAp||F  =  l,n  =  0 
while  ||Ap||f  >  £ 

p"+i  =  (CC'T  +  Xriy^  iCG'^  +  Xrir^-by), 


(a)  (b) 


0  0.2  0.4  0.6  0  0.2  0.4  0.6  0.8 


Fig.  3.  Concentration  estimates  from  simulated  data,  (a)  Bioaerosol, 
(b)  Interferent. 


r-+i  =  (p"+i  +b^-  Ppjxyy/ 
||(p"+i  +  6;f-p,/Ac)+|| 
6;?+!  =6”-bp"+i  -r"+i, 

Ap  =  -  p", 

n  — >■  n  -f  1, 

end 

p  =  p^+^^Dp=  ||/3-poId|lF.Pold  =  P 

end 

Output  p  and  C 


V.  Numerical  Examples 

In  this  section,  we  illustrate  the  multi-aerosol  unmixing 
algorithm  on  simulated  and  actual  release  data  collected  by  the 
U.S.  Army  Frequency  Agile  Lidar  (FAL)  sensor.  The  calcula¬ 
tions  were  done  by  implementing  Algorithm  2  in  MATLAB  on 
a  PC.  The  simulated  FAL  data  were  created  by  injecting  two 
overlapping  aerosol  plumes  having  Gaussian  range  dependence 
both  centered  at  1.5  km  with  width  50  m  into  background  FAL 
data  sets.  Those  sets  consisted  of  transmitted  pulse  waveforms 
and  received  laser  backscatter  from  the  natural  atmosphere  as 
a  function  of  time  step,  wavelength  index,  and  range  cell.  The 
two  injected  plumes  represent  a  bioaerosol  simulant  (bacillus 
BG)  and  interferent  (kaolin  dust).  The  simulant  plume  was 
injected  between  time  steps  200  and  600,  and  the  interferent 
between  time  steps  400  and  800.  A  randomly  varying  peak 
concentration  was  created  from  a  first-order  AR  model.  The 
purpose  of  the  simulator  was  to  generate  data  for  overall  algo¬ 
rithm  performance  verification  including  background  removal 
and  transmitter  pulse  deconvolution.  Because  those  operations 
are  not  relevant  to  the  present  discussion,  we  focus  on  the  com¬ 
parison  of  the  input  concentration  and  backscatter  waveforms 
with  the  estimates  from  the  split  Bregman  algorithm  after  data 
preprocessing. 

Fig.  3  shows  the  concentration  estimates  from  the  unmixing 
algorithm  with  regularization  parameters  set  at  Xr  =  Xc  =  0.05 
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Fig.  4.  Split  Bregman  peak  concentration  estimates  from  simulated  data. 
(Top)  Bioaerosol.  (Bottom)  Interferent. 


Fig.  6.  Concentration  range  dependence  from  simulated  data. 
(Top)  Bioaerosol.  (Bottom)  Interferent. 


Fig.  5.  Bayes  processor  peak  concentration  estimates  from  simulated  data. 
(Top)  Bioaerosol.  (Bottom)  Interferent. 

and  fip  =  He  =  The  plots  show  the  bioaerosol  estimates 
(left)  and  interferent  (right)  versus  time  step  and  range.  The 
calculation  was  initialized  using  the  mean  of  the  backscatter 
spectral  estimates  used  to  train  our  SVM  classifier  for  the 
bioaerosol  and  interferent  classes. 

Fig.  4  compares  the  input  and  estimated  peak  concentration 
over  range  as  a  function  of  time  step  for  the  bioaerosol  (top) 
and  interferent  (bottom).  We  see  that  the  estimates  track  the 
input  throughout  the  run  including  the  plume  overlap  region. 
For  comparison,  Fig.  5  plots  the  peak  concentration  estimates 
from  a  previous  algorithm  based  on  using  backscatter  mean 
spectra  as  priors  in  a  Bayesian  scheme.  Since  the  prior  as  well  as 
data  densities  were  chosen  to  be  Gaussian,  the  Bayes  approach 
in  effect  is  doing  quadratic  norm  processing.  Although  giving 
good  results  when  the  plumes  are  separated  in  time,  this  pro¬ 
cessor  fails  to  produce  good  material  separation  in  the  overlap 
region. 

Figs.  6  and  7  compare  the  concentration  range  dependence 
(true  and  estimated)  and  backscatter  wavelength  dependence 


Fig.  7.  Backscatter  wavelength  dependence  from  simulated  data. 
(Top)  Bioaerosol.  (Bottom)  Interferent. 

at  time  steps  300  and  700  for  the  bioaerosol  and  interferent, 
respectively.  The  results  suggest  that  the  split  Bregman  unmix¬ 
ing  can  provide  both  low  bias  and  low  variance  estimates. 

As  an  example  of  processing  FAL  data  from  a  mixture  of 
aerosols,  we  consider  the  partially  overlapped  release  TD4065 
of  Florida  BG  (a  bioaerosol  simulant)  and  Kuwait  dust  on 
May  31,  2008  in  the  Joint  Ambient  Breeze  Tunnel  (JABT)  at 
Dugway  Proving  Ground,  UT.  The  FAL  sensor  was  located 
about  1 .2  km  from  the  entrance  to  the  JABT.  The  bioaerosol  was 
released  between  time  steps  100  and  476,  and  the  dust  between 
times  261  and  633.  Data  collected  during  the  first  100  time 
steps  were  used  to  estimate  the  backscatter  from  the  natural 
atmosphere.  The  aerosols  were  injected  near  the  far  end  of  the 
tunnel  and  drawn  toward  the  front  of  the  tunnel  by  large  fans, 
where  they  were  exhausted. 

Fig.  8  shows  the  concentration  estimates  from  the  unmixing 
algorithm  on  these  data  with  the  bioaerosol  (dust)  plotted  on 
the  left  (right).  We  see  the  initial  localized  aerosol  appearing  at 
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Fig.  8.  Concentration  estimates  from  JABT  release,  (a)  Bioaerosol.  pjg  iq.  Backscatter  estimates  from  JABT  release,  (a)  Bioaerosol, 
(b)  Interferent.  (b)  Interferent. 
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Fig.  9.  Peak  concentration  estimates  from  JABT  release.  (Top)  Bioaerosol.  Fig-  H-  Classifier  decision  functions  from  JABT  release.  (Top)  Bioaerosol. 
(Bottom)  Interferent.  (Bottom)  Interferent. 


the  indicated  beginning  of  the  releases.  The  subsequent  time 
steps  show  the  plumes  expanding  to  fill  the  space  between 
the  injection  and  extraction  locations.  Fig.  9  plots  the  peak 
concentration  over  range  with  the  bioaerosol  (dust)  shown  at 
the  top  (bottom).  We  see  a  clean  separation  of  the  two  aerosol 
components  in  these  figures. 

The  estimates  of  aerosol  backscatter  are  plotted  in  Fig.  10 
as  a  function  of  time  step  and  wavelength  index  for  the 
16  wavelengths  used  in  these  tests.  Up  to  time  step  100,  we 
see  only  the  random  unit  vectors  produced  by  the  algorithm 
in  the  absence  of  concentration.  After  the  bioaerosol  release 
(shown  on  the  left),  the  estimates  display  a  consistent  structure, 
although  they  are  rather  noisy  between  time  steps  300  and 
450  because  of  low  concentration  levels.  The  corresponding 
dust  backscatter  estimates  (shown  on  the  right)  also  indicate  a 
consistent  structure  throughout  the  dust  release. 

Finally,  the  aerosol  classifier,  a  least-squares  SVM  classifier 
[10],  was  applied  to  the  backscatter  estimates.  Our  classifier 


consists  of  a  three-class  discriminator  for  bioaerosol,  interfer¬ 
ent,  and  null  classes.  The  latter  class  was  included  to  allow  the 
classifier  to  perform  a  detection  function  by  rejecting  spectral 
patterns  that  look  like  neither  bioaerosol  nor  interferent.  Three 
binary  classifiers  were  trained  by  the  1-versus-rest  method 
for  discriminating  samples  from  each  class  against  the  pooled 
alternatives.  At  each  time  step,  the  class  whose  classifier  has  the 
highest  decision  function  value  is  chosen.  Backscatter  estimates 
from  both  material  components  were  processed  this  way,  and 
the  results  for  the  decision  functions  are  plotted  in  Fig.  1 1 .  The 
results  indicate  that  the  classifier  recognized  both  materials  at 
the  correct  times  during  the  partially  overlapped  releases. 

VI.  Summary  and  Conclusion 

The  unmixing  of  elastic  backscatter  data  from  multi¬ 
wavelength,  range-resolved  lidar  for  aerosol  mixtures  has 
defied  solution  by  traditional  regularization  methods  using 
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quadratic  smoothing  constraints.  A  new  shrinkage-based 
Li -regularization  method  for  linear  inverse  problems,  the 
split  Bregman  algorithm,  has  been  successfully  applied  to  the 
aerosol  unmixing  problem.  It  was  illustrated  here  on  simulated 
and  actual  aerosol  mixture  release  data  collected  by  U.S.  Army 
personnel  in  field  testing  at  Dugway  proving  Ground,  UT  using 
the  rapidly  tuned  FAL  sensor. 

The  capability  of  resolving  data  from  mixtures  of  aerosols 
into  their  backscatter  and  concentration  components  is  impor¬ 
tant  to  the  success  of  standoff  sensors  for  biological  detection 
in  realistic  operating  environments  where  interferent  materials 
such  as  smoke  and  dust  will  usually  be  present.  Failure  to 
correctly  perform  this  unmixing  can  lead  to  increased  mis- 
classihcations  that  degrade  the  usefulness  of  potential  sensors 
employing  active  detection. 

Other  possible  standoff  sensing  applications  of  this  unmixing 
technique  include:  1)  the  analysis  of  lidar  data  from  chemical 
releases  where  the  agent  can  have  both  vapor  and  aerosol  phases 
and  interferent  materials  can  include  byproducts  of  an  explosive 
release,  and  2)  the  use  of  thermal  imaging  sensors  operating  in 
the  long  wave  IR  spectra  region. 

Appendix 

Augmented  Lagrangian  Derivation 
OF  Split  Bregman 

Given  the  augmented  Lagrangian  function  L{u,  d,  h)  in  (12), 
we  wish  to  hnd  a  saddle  point  (tt*,  d*,b*)  such  that' 

L{u*,  d*,b)<  L{u*,  d\  b*)  <  L{u,  d,  b*).  (23) 

We  solve  this  saddle  point  problem  by  iterating  between  primal 
and  dual  optimizations 

(it'"”"*'',  =  argminL(M,  d,  (primal) 

u,d 

^  ^k+1  _  ^k+1  (24) 

For  the  primal  minimizations  holding  the  dual  variable  b 
fixed,  we  note  that  because  of  the  introduction  of  the  d  vari¬ 
ables,  the  optimization  over  u  can  be  done  for  fixed  d  by 
differentiation 

dL{u,  d,  b)  _  _  y)  _|_  _l_  _  d'j  (25) 

ou 

to  get 

yk+l  ^  y^k  _  (26) 

Likewise,  although  the  d  variable  optimization  is  still  non¬ 
smooth,  it  is  very  easy  to  solve  by  shrinkage  for  given  u 
and  b 

d^+^  =  S^/^{u'^+^+b^)  (27) 

'This  sketch  follows  the  split  Bregman  treatment  given  hy  Gui-Bo  and 
Xiaohui  [6], 


where  the  shrinkage  function  Sx  is  given  by  (13).  Steps  (26) 
and  (27)  can  be  iterated  for  hxed  b^,  but  we  find  it  more  efficient 
(usually)  to  perform  them  only  once. 
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